Dengue susceptibility analysis

Intro

This analysis aims to show in a naive way how many people are susceptible to dengue and its serotypes in Brazil.

Reading data

Code

pacman::p_load(tidyverse, readxl, geobr)

# Dados do SINAN
denvbr <- read_csv("data/denvbr20152024.csv.gz")

# Populacao CENSO
pop22 <- readxl::read_xlsx(path = "data/POP2022 CD2022_Populacao_Coletada_Imputada_e_Total_Municipio_e_UF.xlsx", 
                    sheet = "Municípios",
                    range = "B3:H5573"
)

rgi <- read_xlsx(path = "data/regioes_geograficas_composicao_por_municipios_2017_20180911.xlsx")

BR0 <- read_municipality()

  |                                                                            
  |                                                                      |   0%
  |                                                                            
  |===                                                                   |   4%
  |                                                                            
  |=====                                                                 |   7%
  |                                                                            
  |========                                                              |  11%
  |                                                                            
  |==========                                                            |  15%
  |                                                                            
  |=============                                                         |  19%
  |                                                                            
  |================                                                      |  22%
  |                                                                            
  |==================                                                    |  26%
  |                                                                            
  |=====================                                                 |  30%
  |                                                                            
  |=======================                                               |  33%
  |                                                                            
  |==========================                                            |  37%
  |                                                                            
  |=============================                                         |  41%
  |                                                                            
  |===============================                                       |  44%
  |                                                                            
  |==================================                                    |  48%
  |                                                                            
  |====================================                                  |  52%
  |                                                                            
  |=======================================                               |  56%
  |                                                                            
  |=========================================                             |  59%
  |                                                                            
  |============================================                          |  63%
  |                                                                            
  |===============================================                       |  67%
  |                                                                            
  |=================================================                     |  70%
  |                                                                            
  |====================================================                  |  74%
  |                                                                            
  |======================================================                |  78%
  |                                                                            
  |=========================================================             |  81%
  |                                                                            
  |============================================================          |  85%
  |                                                                            
  |==============================================================        |  89%
  |                                                                            
  |=================================================================     |  93%
  |                                                                            
  |===================================================================   |  96%
  |                                                                            
  |======================================================================| 100%

Code

data <- pop22 |> transmute(
  UF,
  coduf = `COD. UF`,
  codmun7 = paste0(`COD. UF`,`COD. MUNIC`),
  codmun6 = substr(codmun7,1,6),
  pop = `POP. TOTAL`
)

aux2 <- denvbr |> 
  mutate(SOROTIPO = ifelse(is.na(SOROTIPO),SOROTIPO, paste0("DENV",SOROTIPO))) |> 
  group_by(codmun6 = as.character(ID_MN_RESI)) |> 
  mutate(casos = n()) |> 
  group_by(codmun6, SOROTIPO) |> 
  summarise(
    n = n(),
    casos = casos[1]
  ) |> 
  pivot_wider(values_from = n, names_from = SOROTIPO)
  
data <- data |> left_join(aux2) 

Distribution of the proportion of cases in the last 10 years divided by the 2022 population according to population size.

Code

data |> 
  mutate(
    ratio = casos / pop,
    ratio = replace_na(ratio, 0),
    popsize = case_when(
      pop > 500000 ~ "Pop >500k",
      TRUE ~ "Pop <= 500k"
    )
  ) |> 
  ggplot(aes(x = ratio, color = popsize)) + 
  geom_density() +
  theme_bw()

Serotype

Code

data <- data |> 
  left_join(rgi, by = c("codmun7"="CD_GEOCODI")) |> 
    mutate(
    DENV1 = replace_na(DENV1,0),
    DENV2 = replace_na(DENV2,0),
    DENV3 = replace_na(DENV3,0),
    DENV4 = replace_na(DENV4,0),
    DENVT = DENV1 + DENV2 + DENV3 + DENV4,
    DENVT = ifelse(DENVT==0,1,DENVT)
  ) |> 
  # Calcular proporcoes de sorotipo segundo RGI
  group_by(cod_rgi) |> 
  mutate(
    DENV1p = sum(DENV1) / sum(DENVT), 
    DENV2p = sum(DENV2) / sum(DENVT), 
    DENV3p = sum(DENV3) / sum(DENVT), 
    DENV4p = sum(DENV4) / sum(DENVT) 
  ) |> ungroup() |> 
  # Se a populacao for maior que 100k e ter em 10 anos mais de 100 confirmados
  mutate(
    ratio = casos / pop,
    ratio = replace_na(ratio, 0),
    DENV1p = if_else(pop >= 100000 & DENVT > 100, DENV1/ DENVT, DENV1p), 
    DENV2p = if_else(pop >= 100000 & DENVT > 100, DENV2/ DENVT, DENV2p), 
    DENV3p = if_else(pop >= 100000 & DENVT > 100, DENV3/ DENVT, DENV3p), 
    DENV4p = if_else(pop >= 100000 & DENVT > 100, DENV4/ DENVT, DENV4p)
  ) 

BR <- BR0 |> 
  left_join( data |> 
               mutate(code_muni = as.numeric(codmun7)) |> 
               select(code_muni, ratio, DENV1p:DENV4p)
             )

BR |> 
  mutate(ratio = replace_na(ratio, 0)) |> 
  ggplot(aes(fill = ratio)) + 
  geom_sf(color = NA) + colorspace::scale_fill_continuous_sequential(palette = "inferno") +
  theme_void()

Code

# BR |> 
#   mutate(ratio.cat = cut(ratio, 
#                          breaks = c(0,.2,.5,1.1),
#                          labels = c("<20%", "20-50%", ">50%")
#                          )
#          ) |> 
#   ggplot(aes(fill = ratio.cat)) + 
#   geom_sf(color = NA) + 
  # theme_void()

DENV1

Code

BR |> 
  mutate(ratio = replace_na(ratio, 0) * DENV1p ) |> 
  ggplot(aes(fill = ratio)) + 
  geom_sf(color = NA) + 
  colorspace::scale_fill_continuous_sequential(palette = "inferno") +
  theme_void()

DENV2

Code

BR |> 
  mutate(ratio = replace_na(ratio, 0) * DENV2p ) |> 
  ggplot(aes(fill = ratio)) + 
  geom_sf(color = NA) + 
  colorspace::scale_fill_continuous_sequential(palette = "inferno") +
  theme_void()

DENV3

Code

BR |> 
  mutate(ratio = replace_na(ratio, 0) * DENV3p ) |> 
  ggplot(aes(fill = ratio)) + 
  geom_sf(color = NA) + 
  colorspace::scale_fill_continuous_sequential(palette = "inferno") +
  theme_void()

DENV4

Code

BR |> 
  mutate(ratio = replace_na(ratio, 0) * DENV4p ) |> 
  ggplot(aes(fill = ratio)) + 
  geom_sf(color = NA) + 
  colorspace::scale_fill_continuous_sequential(palette = "inferno") +
  theme_void()

Risk maps

Code

# breaks_risk <- seq(0.5, 1, 0.1)
# labels_risk <- c("0-50%", "60%", "70%", "80%", "90%", "100%")

sf::sf_use_s2(FALSE)

geom_rgi <- rgi |> 
  left_join(BR0 |> 
              mutate(code_muni = as.character(code_muni)), 
            by = c("CD_GEOCODI" = "code_muni")) |> 
  group_by(cod_rgi) |> 
  summarise(geometry = geom |> 
              sf::st_union() |> 
              sfheaders::sf_remove_holes()) |> 
  sf::st_as_sf() |> 
  left_join(rgi)

data_rgi <- geom_rgi |> 
  left_join(data,
            by = c("CD_GEOCODI" = "codmun7"))

denv1.risk <- data_rgi |> 
  mutate(ratio = replace_na(ratio, 0) * DENV1p ) |> 
  ggplot(aes(fill = 1 - ratio)) + 
  geom_sf(color = NA) + 
  colorspace::scale_fill_binned_sequential("inferno")+
  theme_void()+
  labs(title = "DENV1")
denv1.risk

Code

denv2.risk <- data_rgi |> 
  mutate(ratio = replace_na(ratio, 0) * DENV2p ) |> 
  ggplot(aes(fill = 1 - ratio)) + 
  geom_sf(color = NA) + 
  colorspace::scale_fill_binned_sequential("inferno")+
  theme_void()+
  labs(title = "DENV2")
denv2.risk

Code

denv3.risk <- data_rgi |> 
  mutate(ratio = replace_na(ratio, 0) * DENV3p ) |> 
  ggplot(aes(fill = 1 - ratio)) + 
  geom_sf(color = NA) + 
  colorspace::scale_fill_binned_sequential("inferno")+
  theme_void()+
  labs(title = "DENV2")
denv3.risk

Code

denv4.risk <- data_rgi |> 
  mutate(ratio = replace_na(ratio, 0) * DENV4p ) |> 
  ggplot(aes(fill = 1 - ratio)) + 
  geom_sf(color = NA) + 
  colorspace::scale_fill_binned_sequential("inferno")+
  theme_void()+
  labs(title = "DENV4")
denv4.risk

Code

library(patchwork)

patchwork_risk <-(denv1.risk | denv2.risk | denv3.risk | denv4.risk)+
  plot_layout(guides = 'collect')&
  theme(legend.position = "bottom")
patchwork_risk